N91-10 


APRIL 1966 


FLATAU, STEPHENS, DRAINE 


RADIATIVE PROPERTIES OF VISIBLE AND SUBVISIBLB CIRRUS: 
SCATTERING ON HEXAGONAL ICE CRYSTALS 

Piotr J. Flatau, Graeme L. Stephens 

Department of A tmospheric Science , Colorado State University, Fort Collins , CO 808*3 


Bruce T. Drains 

Princeton University Observatory, Peyton Ball, Princeton NJ 08844 
(April 1988) 


1. Introduction 

One of the main objectives of the First In- 
ternational Satellite Cloud Climatology Project 
(ISCCP) Regional Experiment (FIRE) is to pro- 
vide a better understanding of the physics of 
upper level clouds. There are a number of fac- 
tors that have complicated the study of these 
clouds leaving our present understanding of up- 
per level clouds in a less than satisfactory state. 
One factor is simply the distance of the clouds 
from the surface. It is easier to study a low-level 
cloud system because of their proximity to the 
ground. Added to this observational problem is 
the fact that clouds in the upper troposphere 
seem to be governed by an intricate balance of 
not so well understood phenomena. Among the 
factors that have added to the complexity of the 
problem is the fact that particles are no longer 
spherical and both the solid and liquid water 
phases may coexist. Turbulence develops in the 
stable layer environment with its line scales of 
vertical motion. Horizontal eddies (two dimen- 
sional turbulence) may be important for cloud 
morphology, and a host of interactions between 
gravity waves, turbulence, radiation, and micro- 
physics all seem possible. 

This paper concentrates on just one specific 
aspect of cirrus physics, namely on character- 
izing the radiative properties of single, non- 
spherical ice particles. While focused in this 
way, this study provides the basis for further 
more extensive studies of the radiative transfer 
through upper level clouds. Radiation provides 
a potential mechanism for strong feedback be- 
tween the divergence of in-cloud radiative flux 
and the cloud microphysics and ultimately on 
the dynamics of the cloud. 

We will firstly describe some aspects of ice 


cloud microphysics that are relevant to the ra- 
diation calculations. Next, the Discrete Dipole 
Approximation (DDA) is introduced and some 
new results of scattering by irregular crystals 
are presented. We also adopt the Anomalous 
Diffraction Theory (ADT) to investigate the 
scattering properties of even larger crystals. In 
this way we are able to determine the scatter- 
ing properties of non spherical particles over a 
range of particle sizes. The study reported here 
is still preliminary and at the time of writing 
this abstract, the results are incomplete. We 
aim to incorporate the microphysics data col- 
lected during FIRE into these calculations and 
in a companion study plan to use these scat- 
tering properties to determine cloud radiative 
heating rates. 

2. Cloud microphysics 

Characterizing the shape and size of ice crys- 
tals in terms of their environment continues to 
be a subject of extensive research. For envi- 
ronments typical of cirrus clouds and for even 
the colder environments of polar stratospheric 
clouds (PSC’s), ice exists in single- and poly- 
crystalline forms. The more commonly per- 
ceived large composite crystals, of the order of 
500/zm, are actually observed more readily at 
lower temperatures. Crystal sizes of the order 
of 10 fim and Ifim are, respectively, more repre- 
sentative of high cirrus in the tropics (Heyms- 
field, 1986) and PSC’s (Rosen e* ol, 1988). To 
emphasize this point we present the variation of 
crystal shape drawn in proportion to their size 
in Fig. 1 over the temperature range — 20°C to 
— 90°C. The crystals at cirrus temperatures, say 
— 40°C, are hexagonal columns although varia- 
tions of this type of crystal exist in the form of 
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hexagonal bullets and hollow columns. In this 
study we will focus only on single solid hexag- 
onal crystals although the foundation to study 
crystals of more complex shape has been devel- 
oped. 



Fig. 1 Typical sizes of ice crystals encountered 
in the temperature range (-20, — 90) °C. The 
crystal sizes are approximately drawn to scale . 
Notice that the warm cirrus particle (all hatched 
in the low part of the diagram) is so large that 
only part is presented. 


3. Integral formulation of Maxwell equations 
and its approximation 

Only a very brief outline of the scattering 
methods are described here. We begin with 
Maxwell equations written in the integral form 
(Saxon, 1955) 

E(x) = E 0 (x) + 4nk 2 

X / G(x,x') • ~ 1 • E(x') d V (1) 

J 4jr 

where Eo(x) is the incident electric field, e(x) 


is a dielectric tensor, and Green’s function is 
defined as 


a ikr 


G(x,x') = (l-rr) — 


+ (3rr 1) ( fc2r2 jbj .) 


gtfcr 

47r r 


( 2 ) 


where r = x — x? t r = |x — x'|, and r = r/r. 
Boldfaced 1 is a unit matrix, and rr is a dyadic 
(or tensor) multiplication of two unit vectors r. 
The incident electric field Eo(x) is 

E 0 (*) = E 0 ee iki * (3) 

where e is a unit vector specifying the wave’s 
polarization, E 0 is the intensity of incident 
wave, k is a unit vector defining the wave 
propagation, k == 2tt/A, and A is wavelength. 
We now seek a solution of this integral equa- 
tion to obtain the electric field and all the 
relevant scattering properties follow from this 
solution. Unfortunately, the integral (Fred- 
holm) equation (1) is not easily solvable al- 
though it is amenable to a variety of approx- 
imations commonly used in scattering prob- 
lems. In fact the approximate methods such 
as the Rayleigh approximation (also known as 
the Rayleigh-Gans- Debye (RGD), the Rayleigh- 
Gans-Rocard or the Kirchoff or the first Bom 
approximation), the second and higher Bom 
approximations, Anomalous Diffraction Theory 
(ADT) also known as High Energy Approxima- 
tion (HEA) or the equivalent WKB method, 
the Discrete Dipole Approximation (DDA) and 
the Extended Boundary Method can all be de- 
rived directly from (l) (Saxon, 1955; Flatau and 
Stephens, 1988b). 


4. Discrete Dipole Approximation 

A useful and general way to solve (1) is to dis- 
cretize the integral term thus reducing the in- 
tegral equation to a linear system of equations. 
In the DDA approximation one assumes that 
each volume element may be replaced by a point 
dipole, whose (tensor) polarizability is related 
to the volume of the element and the dielectric 
tensor using the Clausius-Mossotti relations, 
corrected for the effects of radiative reaction. 
Each of the dipoles acquires a (time- dependent) 
polarization in response to the electric field at 
the position of the dipole, which includes con- 
tributions from all of the other dipoles plus the 
incident wave. This is the approach of the DDA 
as originally formulated by deVoe (1964) and 
Purcell and Pennyp acker (1973). Consider our 
hexagonal crystal in this instance divided into 
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N discrete domains of volume Vj (Fig. 2) cen- 
tered around point xj, where j = 1 
The volume of such sin elementary domain is 


Vj = 


V 

N 


( 4 ) 


where V is the volume of the particle, and N 
is the number of domains. Thus we are able to 
form a set of 3 N equations from (1) since there 
are three components of the B field at each x*. 
The resulting equation set can then be repre- 
sented in the form 


AV = € 0 (5) 


where P is the 3N- dimensional vector whose el- 
ements are the dipole moments of the N dipoles. 
For further details of the matrix structure and 
computational solution procedures the inter- 
ested reader can consult the recent work of 
Draine (1988). The electric field everywhere in 
the three-dimensional space is given in terms 
of the structure matrix A which depends only 
on particle’s shape, refractive index and inci- 
dent wavelength but is independent of the direc- 
tion and polarization of the incident wave Eo(x). 
Once V is determined then single scattering 
properties of the particle can be derived from 
the relationships reported by Draine (1988). 



Fig. 2 Hexagonal discrete dipole arrays with 
N — 384. Higher resolution was actually used, 
see text for further comments . 

The DDA method has been applied to calcu- 
late the scattering properties of small hexago- 
nal ice crystals. Figure 2 presented the discrete 
dipole array for N — 384 dipoles but for the cal- 
culations reported in this paper actually used 
2208 dipoles. The wavelength of X = 3.8pm 
and the corresponding value for the refractive 
index of ice m = 1.38 + *0.0067 was employed 


as was A == 10.8 pm with a corresponding value 
m = 1.089 + *0.182. The wavelengths rep- 
resent near infrared and infrared regions with 
small and relatively large absorption, respec- 
tively and correspond to central wavelengths of 
the AVHRR imager. The 10.8 pm wavelength 
is also a relevant wavelength for application to 
CO 2 lidar studies. 

The scattered intensity is presented in Fig. 3 
as a function of scattering angle for two scat- 
tering planes, one perpendicular to E 0 and the 
other parallel to Eq. 



Fig. 3 Intensities for hexagonal ice for an in- 
cident wave perpendicular to the crystal (aver- 
age) and A = 3.8 pm. Results for two scattering 
planes are presented . Crystal length L = 8pm, 
crystal radius a = 2 pm. N = 2208 discrete 
dipoles were used . 

The incident direction was taken to be normal 
to the L-axis of the hexagon, the length of the 
crystal used in the calculations is specified by 
L % and the aspect ratio is defined as p = L/2a 
where 2a is the width of the crystal. Results are 
shown for two cases: incident electric field par- 
allel to, and perpendicular to, the L-axis of the 
hexagon. Because the hexagon is not rotation- 
ally symmetric, the scattering problem depends 
upon the orientation of the hexagon. Two cases 
have been considered: (1) with the hexagon ori- 
ented so that the k vector is normal to one of the 
6 rectangular faces, and (2) with the hexagon 
rotated by 30 degrees, so that the k vector is 
parallel to a line connecting a vertex to the cen- 
ter of the hexagon. (Each of these two cases 
has reflection symmetry wliich greatly reduces 
the amount of computing required to obtain a 
solution). The scattering intensities shown in 
Figures 3 and 4 are the averages of the scatter- 
ing intensities for the above two cases. 
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Fig. 4 Same as Fig . 3 but for A = 10.8/zm. 
Crystal length L — 8/im, crystal radius a = 
2/xm. 

Figure 5 shows the results of the Mie calcula- 
tion for a sphere with a volume equal to that of 
the hexagonal crystal. The intensities |/|jj and 

j/l^ are presented here as a function of scatter- 
ing angle and 


i/i ! = f 


dC,ca 

da 


( 6 ) 


(see e.g. Bohren said Huffman, 1983). There sire 
several features worthy of mention in comparing 
between Figs. 3, 4 and 5. For example the 
backscatter linear polarization, defined as 


6(6 


180) = l/lK g = 180) - |/|jj(fl = 180) 
|/li(0 = 180) + |/|j(0 = 180) 


( ? ) 

is non- vanishing in case of hexagonals as com- 
pared to a zero depolarization for spheres. This 
is a well known characteristic of scattering by 
non-spherical particles but the quantitative re- 
lationship between 6 and particle shape has 
never been convincingly established. The first 
minimum in intensity doesn’t correspond to the 
Mie case, thus indicating that care has to be 
taken when interpreting results from the for- 
ward scattering PMS probes. This again is well 
established experimentally fact but its theoret- 
ical confirmation has to date been scarce. It 
seems that the extinction can be modelled, as 


it usually is, using a sphere of equivalent vol- 
ume or surface, and we plan to perform more 
detailed calculations for the angularly averaged 
case. 



Fig. 5 Mie results for equivalent volume spheres 
r = 3.41/im, and A = 3.8 jun. 


5. Anomalous Diffraction Theory 

The anomalous diffraction theory (ADT) be- 
longs to the broad class of approximations (to- 
gether with the Rayleigh-Gans-Debye or Bom 
expansions for example) in which one makes 
some assumption about the electric field under 
the integral of (1). This was noticed by Saxon 
(1955) and derived independently by van de 
Hulst (1957) on the basis of geometrical optics 
and diffraction. In the present context one as- 
sumes that the electric field inside the particle 
is given by 


E(x) = Eoeexp [tfck • (b + e z z) - x(b,z)j (8) 

where x(b, z) corresponds to the phase change 
due to the cnanged refractive index inside the 
particle. The impact vector b is perpendicu- 
lar to the z axis (see Fig. 6). The explicit 
formulation of x(b, z) and relation of ADT to 
other theories such as High Energy Approxi- 
mation and DDA is contained in Flatau and 
Stephens (1988). The anomalous diffraction 
theory (ADT) (Stephens, 1984) holds for par- 
ticles with a refractive index close to unity 
(m — 1) 1 and for a particle with a size to 
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wavelength ratio z > 1< Since the index of re* 
fraction is close to 1, the problem of total extinc- 
tion in the ADT theory reduces to calculations 
of the interference between the almost straight 
transmission and the light diffracted according 
to Huygens’ principle. 



Fig. 6 Geometry for Anomalous Diffraction 
Theory calculations, b is an impact parameter, 
t the position vector, k the direction of plane 
wave. The shadow area is in the far field and 
only forward scattering is considered here. 



Fig. 7 Extinction efficiency (convening or os- 
cillating around 2 for large phase shift values) 
and absorption efficiency for hexagonal ice crys- 
tals for 3 different values of complex refractive 
index as a function of phase shift. 

Thus, in the anomalous diffraction approxi- 
mation, the forward scattering amplitude S(0) 


is given by 

*<*-*>- £/>-'■*■)" (9) 
where — kd(m — 1), d is the particle thick- 
ness, and m is the complex refraction: m = 
n - in'. The quantity <f>* is the complex phase 
shift of light passing through the particle rela- 
tive to that passing around it. The kdn' term 
contributes to absorptive attenuation. 

The extinction ana absorption coefficients are 
defined as 

C eet e"** ) dA (10) 

and 

Caht = ( 1_ e ‘ iWn ) dA < u ) 
and the single scattering albedo is 

Oext ~ 

" “ Cext % 

Therefore application of the ADT requires the 
relatively straight forward determination of the 
path length through the particle and evaluation 
of integrals for C** and C^. Notice that (8) is 
more general because it provides, in principle, 
the full phase function. Calculations of C exty 
Cabs, and u> for hexagonal crystals are presented 
in Figs. 7 and 8 using this approximation for the 
geometry depicted in Fig. 6. Full details of the 
method are planned in a forth coming paper and 
we also plan to compare these results, including 
phase functions, with solutions obtained from 
geometric optics and DDA. 




Fig. 8 Single scattering albedo for hexagonal 
ice crystal in ADT approximation for 3 different 
values of complex refractive index . 
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6. Summary 

This paper reports on new calculations of the 
scattering of radiation by hexagonal ice crys- 
tals. Results are presented for small crystals 
which are not only the first of their type but also 
employ methods, to our knowledge, not previ- 
ously used in atmospheric scattering problems. 
Scattering properties of large hexagonal crystals 
were also modelled using more approximate the- 
ories. While the results presented in this paper 
apply only to hexagonal crystals, the methods 
are general and particles of other geometries will 
be considered in further studies. 

The work reported in this paper is prelimi- 
nary and we plan to use the microphysic data 
collected during FIRE into these calculations 
to obtain the optical properties needed in the 
radiative transfer simulations of the cloud flux 
measurements. 
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